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We study single hole dynamics in the bilayer Heisenberg and Kondo Necklace models. Those 
models exhibit a magnetic order-disorder quantum phase transition as a function of the interlayer 
coupling J±. At strong coupling in the disordered phase, both models have a single-hole dispersion 
relation with band maximum at p — (n, n) and an effective mass at this p— point which scales as 
the hopping matrix element t. In the Kondo Necklace model, we show that the effective mass at 
p — (tt, 7r) remains finite for all considered values of J± such that the strong coupling features of the 
dispersion relation are apparent down to weak coupling. In contrast, in the bilayer Heisenberg model, 
the effective mass diverges at a finite value of J± . This divergence of the effective mass is unrelated 
to the magnetic quantum phase transition and at weak coupling the dispersion relation maps onto 
that of a single hole doped in a planar antiferromagnet with band maximum at p = (tt/2, tt/2). We 
equally study the behavior of the quasiparticle residue in the vicinity of the magnetic quantum phase 
transition both for a mobile and static hole. In contrast to analytical approaches, our numerical 
results do not unambiguously support the fact that the quasiparticle residue of the static hole 
vanishes in the vicinity of the critical point. The above results are obtained with a generalized 
version of the loop algorithm to include single hole dynamics on lattice sizes up to 20 x 20. 
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I. INTRODUCTION 

The modeling of heavy fermion systems is based on an 
array of localized spin degrees of freedom coupled anti- 
ferromagnetically to conduction electrons. Those mod- 
els show competing interactions which lead to magnetic 
quantum phase transitions as a function of the antifer- 
romagnetic exchange interaction J. Kondo screening of 
the localized spins, dominant at large J, favors a param- 
agnetic heavy fermion ground state, where the localized 
spins participate in the Luttinger volume. In contrast, 
the RKKY interaction favors magnetic ordering and is 
dominant at small values of J. There has recently been 
renewed interest concerning the understanding this quan- 
tum phase transition. In particular, recent Hall experi- 
ments pj suggest the interpretation that in the vicinity 
of the quantum phase transition the localized spins drop 
out of the Luttinger volume. Starting from the param- 
agnetic phase, this transition from a large to small Fermi 
surface should coincide with a effective mass divergence 
of the heavy fermion band. 

Motivated by the above, we consider here a very simpli- 
fied situation namely that of a doped hole in the Kondo 
insulating state as realized by the Kondo necklace and 
related models. Although this is not of direct relevance 
for the study of the Fermi surface, it does allow us to in- 
vestigate the form of the quasiparticle dispersion relation 
from strong to weak coupling for a variety of models. Our 
aim here is two fold. On one hand we address the ques- 
tion of the divergence of the effective mass as a function 
of coupling for different models, and on the other hand 
the fate of the quasiparticle residue in the vicinity of the 
quantum phase transition. 

The KLM emerges from the periodic Anderson model 



(PAM), where we have localized orbitals (LO) with on- 
site Hubbard interaction Uf and extended orbitals (EO), 
which form a conduction band with dispersion e(p) = 
—2t (cosp x + cospy). The overlap between the LOs and 
the EOs within each unit cell is described by the hy- 
bridization matrix element V. For large Uf charge fluc- 
tuations on the localized orbitals becomes negligible and 
the PAM maps via the Schrieffer- Wolff transformation 
onto the KLM [10: 
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Here 5j and are spin 1/2 operators for the extended 
orbitals and the localized orbitals respectively. In the 
first term, which represents the hopping processes, the 
fermionic operators £p a (c pcr ) create (annihilate) elec- 
trons in the conduction band with wave vector p and 
z-component of spin a. At half-filling - one conduc- 
tion electron per localized spin - the two-dimensional 
KLM is an insulator and shows a magnetic order-disorder 
quantum phase transition at a critical value of J c jt = 
1.45 ±0.05 0. 

By taking into account an additional Coulomb repulsion 
U between electrons within the conduction band, one ob- 
tains a modification of the KLM, the [/KLM: 
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Here, fn a = c] CT Ci CT is the density operator for electrons 
with spin a in the conduction band. The additional 
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FIG. 1: (a) Isotropic bilayer Heisenberg model, (b) Kondo 
Necklace model, that is related to the [/KLM. In both cases 
the system dimerizes for large J±, so that the AF ordering 
breaks down. 

Coulomb repulsion displaces the quantum critical point 
towards smaller value of J c /t. However the physics, in 
particular the single hole dynamics, remains unchange 
[5j- This allows us to take the limit U/t — > oo to map 
the UKLM onto a Kondo necklace model (KNM) which 
we write as: 

H = JiEW + EE^T- ( 3 ) 
» (»i> ™ 

»(m) 

Here is a spin 1/2 operator, which acts on a spin 

degree of freedom at site i. stands for the intralayer 
exchange and the upper index m — 1,2 labels the two 
different layers. The interlayer exchange, formerly the 
AF coupling J between LOs and EOs, is now character- 
ized by Jj_. Clearly, since we have motivated the KNM 
from a strong coupling limit of the UKLM, we have to 
set: 

ee J,, J* 2) = for the KNM. (4) 

The above models all have in common that the only inter- 
action between the localized spins stems from the RKKY 
interaction. This in turn leads to the fact that at J = 
for the KLM and UKLM or J± = for the KNM the 
ground state is macroscopically degenerate. To lift the 
pathology we finally consider a Bilayer Heisenberg Model 
(BHM) in which an independent exchange between the 
localized spins is explicitly included in the Hamiltonian. 
Hence we will equally consider an Isotropic BHM which 
takes the form of Eq. 10) with: 

= J,{ 2) ee J|| for the isotropic BHM. (5) 

Both the KNM and BHM systems are sketched in FIG. 

m 



The main results and organization of the paper are the 
following. In section ^ we give a short overview of the 
quantum Monte Carlo (QMC) method. We use a gener- 
alization of the loop algorithm which allows for the cal- 
culation of the imaginary time Green's function of the 
doped hole |6|. Dynamical information is obtained with 
a stochastic Maximum Entropy method 0, @- I n the 
first part of section IIIII we present our results for the 
spin dynamics. This includes the determination of the 
quantum critical point for the isotropic BHM as well as 
the Kondo Necklace model (KNM) by QMC methods. 
In the second part of that section we analyze the sin- 
gle particle spectral function. It turns out, that there 
are significant differences between the models. We can 
identify two classes of models: In the isotropic BHM the 
dispersion is continously deformed with decreasing inter- 
planar coupling J±/ J|| resulting in a displacement of the 
maximum fromp = (7r,7r) top = (£, ^). In other words, 
the effective mass - as defined by the inverse curvature 
of the quasiparticle dispersion relation - at p = (tt, tt) di- 
verges at a finite value of the interplanar coupling. This 
divergence of the effective mass is not related to the mag- 
netic order-disorder transition. In contrast, in the KLM 
related models, UKLM and KNM, the maximum of the 
quasiparticle dispersion relation is pinned at p = (tt, tt) 
irrespective of the value of the interplanar coupling. In 
those models the effective mass at p = (tt, tt) grows as a 
function of decreasing interplanar coupling, but remains 
finite. 

In section Hvl we turn to the analysis of the quasi particle 
residue (QPR) across the quantum phase transition. To 
gain intuition, we first carry out an approximate calcula- 
tion in the lines of Ref . . The physics of the spin sys- 
tem may be solved in the framework of a bond mean-field 
calculation. Here, the disordered phase is described in 
terms of a condensate of singlets between the planes and 
gaped spin 1 excitations (magnons). At the critical point 
the magnons condense at the AF wave vector thus gen- 
crating the static antiferromagnetic order. Within this 
framework one can compute the coupling of the mobile 
hole with the magnetic fluctuations and study the hole 
dynamics within a self-consistent Born approximation. 
The result of the calculation shows that the quasiparti- 
cle weight at wave vectors on the magnetic Brillouin zone 
[e(p) = e(p + Q) with Q = (tt, tt) ] vanish as the square 
root of the spin gap. In contrast the QMC determina- 
tion of the quasiparticle residue on lattices up to 20 x 20 
for static and dynamical holes does not unambiguously 
support this point of view. 



II. NUMERICAL METHODS 

We use the world line QMC method with loop updates 
[TrU to investigate the physics of the BHM and KNM. To 
investigate the spin dynamics we compute both the spin 
stiffness as well as the dynamical spin structure factor. 
Our analysis of the single hole dynamics is based on the 
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calculation of the imaginary time Green's function. Ana- 
lytical continuation with the use of the stochastic Maxent 
Method provides the spectral function and the quasipar- 
ticle residue is extracted from the asymptotic behavior of 
the imaginary time Green's function. Below, we discuss 
in more details the calculation of each observables. 



Spin Stiffness 

To probe for long-ranged magnetic order we introduce a 
continuous twist in spin space which, when cumulated 
along the length L along (e.g.) the cc-axis, amounts to a 
twist of angle <j> around a certain spin axis e. This means 
thus the boundary conditions read: Si+Le x = R [e, 4>] Si, 
where R [e, 4>] is a matrix describing an SO(3) rotation 
around the axis e by the angle (j>. The spin stiffness is 
then defined as 



1 



Ps = 



L 



1 d 2 



(0) 



0=0 



with (3 as inverse temperature, L as the linear size of the 
system, d the dimensionality and Z((f>) the twist depen- 
dent partition function. In the presence of long-range 
order p s takes a finite value and in a disordered phase it 
vanishes. 

Within the world-line algorithm, the spin stiffness is re- 
lated to the winding number W x of the world line con- 
figurations along the axis of cumulatively twisted spins 
(e.g. x-axis). In particular, in the limit At — > it takes 
the simple form 



L d (3 x 



Spin Correlations 



(7) 



Within the QMC it is easy to obtain the spin correla- 
tions (S-(T)Sj(0)) in real space and imaginary time r, 
where the imaginary time evolution of the spin operator 



reads S% (r) = e 



tH 



S*e 



-tH 



Its representation in momen- 



tum space is related to the dynamical spin susceptibility 
S(q, ui) via: 



(S z q {r) 



S z - q (0)) = - 

7T 



dwe-™S(g,a;). 



(8) 
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tiy using the Stochastic Maximum Entropy (ME) method 
we can extract the dynamical spin susceptibility. For 
large r the spin correlation function is dominated by the 
lowest excitation: 



lim {S:{t)S'(0)) cx e 



-n(q)r 



(9) 



where £l(q) stands for momentum dependent gap to the 
first spin excitation. Thus, we obtain the gap energy A 
from the asymptotic behavior of the spin correlations: 
A = min [0(g)] . 



The Green's Function 

To incorporate the dynamics of a single hole into the 
KNM and BHM, we consider the t J-model 

Ht.j = V s [- +c) .c i0 -) 

which describes the more general case of arbitrary filling. 
Here, i and j denote lattice sites of the bilayer BHM, tij 
the hopping amplitude, Jij the exchange, fij = cl a Ci a , 
and the sums run over nearest inter- and intraplane 
neighbors. Finally Vs is a projection operator onto the 
subspace S with no double occupation. We apply a map- 
ping, introduced by Angelucci [111], which separates the 
spin degree of freedom and the occupation number. 
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(ii) 



fj and fi are spinless fermion operators which act on the 
charge degree of freedom and create (annihilate) a hole 
at site i: fhl, c) = |0, a), af are ladder operators for the 
spin degree of freedom and <r?' = i(l±<7?) are projector 
operators acting on the spin degree of freedom. Within 
this base the Hamilton of the t J-model (|1(J|) writes: 



H tJ = V s [ [fjfiPij + h.c] 



2 "^/(Ai l)^ij 



Vs (12) 



where Pi. 



\{0i0j 



ij ~ 2^ uiu 3 1) anC l ~ 1 — fi fi 

Vs = Yli (l — fi fi^i^t) 1S a projection operator in An- 
gelucci representation which projects into the subspace 
S. This representation (|12fl has two important advan- 
tages which facilitate numerical simulations: (i) Because 
the Hamiltonian commutes with the projection operator: 
[H t j,Vs] — 0, the bare Hamiltonian (H t j without pro- 
jections) generates only states of subspace S provided 
that the initial state is in the relevant subspace. (ii) The 
Hamiltonian is bilinear in the spinless fermion operators. 
Within the Angelucci representation the Green's function 
reads: 



Gji(r) = (^• + (r)/ i (r)af' + (0)//(0)) 



(13) 



The time evolution 



'(r)fj(r) 



imaginary time is given by: 
'fje- T ^ J . The authors of Ref. 



show in details how to implement the Green's func- 
tion into the world line algorithm of our QMC simula- 
tion. The spin dynamics is simulated with the loop algo- 
rithm. For each fixed spin configuration, one can readily 
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compute the Green's function since the Hamiltonian is 
bilinear in the spinless fermion operators /. 
From the Green's function G p (r) we can extract the sin- 
gle particle spectral function A(p, uj) with the Stochastic 
Maximum Entropy: 



GJr) = - [ dLoe- TLU A(p,-uj) 
n Jo 



(14) 



In the T = limit the asymptotic form of the Green's 
function reads: 



G p (r) = \(^- 1 \c p \^)re 



(15) 



where /i is the chemical potential. As apparent, the pref- 
actor, 



z P = \m \cp\n 



i.N\\2 



(16) 



is nothing but the quasiparticle residue. Hence from the 
asymptotic form of the single particle Green's function, 
we can read off the quasiparticle residue. 



III. SPIN AND HOLE DYNAMICS 

In this section we present our results for the spin dynam- 
ics as well as for the spectral function of a doped mobile 
hole. 
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FIG. 2: (a) Spin stiffness p 3 as a function of linear lattice size 
L for different interplanar couplings J±/J\\ in the KN model. 
Extrapolation to the thermodynamic limit is carried out by 
fitting to the form a + b/L (b) Extrapolated value of the spin 
stiffness as a function of J±/Ju. The dashes line corresponds 
to the fit according to the form of Eq. QZp . 



A. Spin Dynamics 

All considered models, KLM, UKLM, KNM and BHM, 
show a quantum phase transition between an antiferro- 
magnetic ordered phase and a disordered phase. It is 
believed, that all models belong to the same universality 
class. To demonstrate this generic property and to test 
our numerical method we determine the quantum critical 
point as well as critical exponents in the isotropic BHM 
and KNM. Fig. plots the spin stiffness for the KNM 
as a function of lattice size. The extrapolated data is 
plotted in Fig. 03. We fit the data to the form: 



(17) 



to obtain (J±/Jn) c — 1.360 ± 0.017 and a critical expo- 
nent of v — 0.582 ± 0.077, whic h ag rees (within the error 
bars) with the value of Ref. v = 0.685 ± 0.035. 

Similar data for the BHM localizes the quantum criti- 
cal point at (J±/J\\)c = 2.5121 ± 0.0044, which conforms 
roughly the literature value (Jx/J^f = 2.525 ± 0.002 
of Ref. 01 . For the critical exponent we obtain v = 
0.7357 ± 0.044. Again this is in good agreement with 
the critical exponent specified in Refs. ^2[- I* 1 -f^f- E-l 
the BHM and the KNM are observed by dimer series ex- 
pansions. Within this framework our numerical results 
are reflected quite well. FIG. |21 plot the dynamical spin 



structure factor as a function of J±/J\\ for the BHM. In 
the deeply disordered phase the dispersion has a cosine- 
like shape. In the limit Jj_ — ► oo the ground state wave 
function is a tensor product of singlets in each unit cell. 
Starting from this state, a magnon corresponds to break- 
ing a singlet to form a triplet. In first order perturbation 
theory in J± / J\\ , the magnon acquires a dispersion rela- 
tion: 

n( 9 ) « Jl + \ J|| 7 (g) (18) 

with 7(g) = 2 {cos(q x ) + cos(q y )). This approximative 
approach is roughly consistent with the large- J± case in 
Fig. Et- As as function of decreasing coupling J± the spin 
gap progressively closes (see Fig. and at the critical 
coupling the magnons at q = (ir, ir) condense to form the 
antifcrromagnetic order. This physics is captured by the 
bond mean field approximation which we discuss below. 



Bond Operator Mean Field Approach 

The bond mean field approach ^| is a strong coupling 
approximation in J±. The spins between layers domi- 
nantly form singlets and the density of triplets is "low". 
This assumption allows one to neglect triplet-triplet in- 
teraction. The bond operator representation describes 
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(a) JJJ t = 3.5 
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FIG. 3: Dynamical spin susceptibility, respectively magnon 
dispersion for different coupling ratios on a 12 x 12 square 
lattice. 



the system in a base of pairs of coupled spins, which can 
either be in a singlet or triplet state. 

1 



\a)i = §l\0)i = -j={\ U)i ~ I it)*) 

i*x>* = 4io)i = ^(i m<-iu>o 



\ty)i = tl v \0)i 
\t*)i 



iy 1 



i 

i 



(I TT)< + I U><) 



4l0)i = ^(l U)i + 1 IT)») 



(19) 



The operators tr and satisfy Bose commutation rules 
provided that we impose the constraint 

4*+5>l*.-a = 1 ■ ( 2 0) 

ot 

Since the original spin 1/2 degrees of freedom reads, 

&ia "* = \{^\ha ± — Cap-ftlpiiy) , (21) 

/3 7 



the Hamiltonian can be rewritten in the bond opera- 
tor representation as: 



i 



IV? I ) 



<*i> Q 

a,/3,7 

is a Lagrange parameter which enforces locally the 
constraint (|20|> . The interplanar part shows the charac- 
teristic Hamiltonian of two antiferromagnetically coupled 
spins whereas the intraplanar part includes the interac- 
tion between singlets and triplets of different bonds. We 
now follow the standard method of Sachdev and Bhatt 
[l5|. In the disordered phase we expect a singlet con- 
densate (s = (s) 0) and impose the constraint only on 
average (pn = fi) . As mentioned above we neglect triplet- 
triplet interactions. Apart from a constant we obtain the 
following mean field Hamiltonian in momentum space: 



Hmfa - EE^'V?" 



a q 



EE 2 ^°^-i a 



h.c. 



where 



A„ 



J± 
4 



B q = J\\s 



- H + J\\s 2 (cos(g x ) + cos(q y )) 
2t cos(q x ) + cos(q y )) . 



(22) 

(23) 
(24) 



The parameter /i and s — (s) are determined by 
the saddle-point equations: {OHmfa/ dn) = and 
(dHMFA/ds) = 0. The Hamiltonian is diagonalized by 



a Bogoliubov transformation: a qa 



V q t-q a . Ill 



terms of magnon creation and annihilation operators the 
Mean field Hamiltonian IL'l't writes: 



Hmfa = E E ^l^l^i 



(25) 



The Bogoliubov coefficients u q and v q satisfy the relation 
ui — Vg = 1, which follows from the bosonic nature of the 



magnons: [a q ,a q , 



The coefficients are given by 



1 q> v q — 



I A q 1 
2Q(q) 2 ' 



(26) 



where £l(q) = J A q — B q is the magnon dispersion. In 
the vicinity of the critical point it can be approximated 
by 



(27) 
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FIG. 4: (a) Spin correlation function (J_i_/J\\ = 2.7) at p = 
(tt, tt) for a 12 x 12 lattice in the BHM. (Inverse temperature 
/3 J|| = 30.0, At J|| = 0.02) (b) Spin gap A at p = (tt, tt) for 
different coupling ratios Jj_/Jii. The data for a 12 x 12 lattice 
is fitted by A oc (g — 9c)~ z " with g = J±JJ\\ and literature 
values p c = 2.525 ± 0.002 and z = 1 



with A the energy gap to magnon excitations, v s the 
magnon velocity and Q = (tt, tt). Eq. I|27l) gives an accu- 
rate description of the dispersion relation in the vicinity 
of the critical point (see FIG. |3J:). At the critical point 
the gap A vanishes, so that the triplets can condense thus 
forming the AF static ordering. 



B. Hole dynamics 

We now dope our systems with a single mobile hole and 
restrict its motion to one layer thereby staying in the 
spirit of Kondo lattice models. To understand the cou- 
pling of the hole to magnetic fluctuations within the mag- 
netic disordered phase we can extend the previously de- 
scribed bond mean-held approximation (See Eq. J2HJ)) to 
account for the hole motion. For this we introduce the 
operator h, (hi a ), that creates (anhilates) a hole with 
spin a in layer 1 at site i. 

hl\vac) = \0a)i (28) 

|ciC2)i denotes a dimer state at site i with spin <j\ in 
layer 1 and spin <ji in layer 2. The Hamiltonian now 



writes [Til ]: 

H = Y,Sl(q)ala q + Y,e(p)hlh p (29) 

q p 

+ 9& q ^ a 9 ■ (K+q ah p) + h - c - 
pq 

with spinor h p — {h p ^,h p {) and vector 6t q = 
). a = (a 1 , a 2 , a 3 ) denotes the Pauli ma- 
trices. The coupling strength between the hole and 
magnons is given by g(p 7 q). We discuss g(p,q) in de- 
tail later in section llVl For the bare hole dispersion the 
calculation yields 

e(p) = +ts 2 (cos(p x ) + cos(p y )) . (30) 

In the limit Jj_ — > oo the magnon excitation energy 
diverges (see Eq. and hence the coupling of the 

hole to magnetic excitations becomes negligible. In this 
limit the magnon excitations become quite rare, so that: 
,s = (s) ~ 1. Thus, in the strong coupling region we 
obtain from (|29|) a hole dispersion relation: 

E(p) = t(cos( Px ) + cos(p. y )) . (31) 

This agrees with the result given by applying perturba- 
tion theory in t/J± 0. 

As apparent from Figs. [S] and this strong coupling 
behavior is reproduced by the Monte Carlo simulations 
where the dispersion exhibits a cosine form with maxi- 
mum at p — (tt, tt) . The form of this dispersion relation 
directly reflects the singlet formation - in other words 
Kondo screening - between spin degrees of freedom on 
different layers. We note that this strong coupling be- 
havior of the dispersion relation sets in at larger values 
of Ji/J|| for the BHM than for the KNM. This is quite 
reasonable since in the BHM the single bonds are coupled 
among each other within both layers. 
With decreasing coupling ratio the bandwidth of the 
quasiparticle dispersion relation diminishes but the over- 
all features of the strong coupling remain. 
In the weak coupling limit we observe considerable differ- 
ences between the single particle spectrum of the BHM 
and KNM. Let us start with the BHM. For this model 
the point J±/J\\ — is well defined (i.e. the ground is 
non-degenerate on any finite lattice) and corresponds to 
two independent Heisenberg planes with mobile hole in 
the upper plane. The problem of the single hole in a two 
dimensional Heisenberg model has been addressed in the 
framework of the self-consistent Born approximation > 
and yields a dispersion relation given by: 

E(p) = J,, (cos(p)s + cos( Py )f . (32) 

Since at J±/J\\ = we have a well defined ground state 
we can expect that turning on a small value of J±/J\\ 
will not alter the single hole dispersion relation. This 
point of view is confirmed in Fig. |SJ At J±/J\\ = 1, 
the single hole dispersion relation follows of Eq. I|32|l . 
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(a) J ± /J l{ = 10.0 
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(b) Jj/Jj, = 2.4 
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(c) Jx/Jy = 2.0 
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(d) Jx/Jy = 1.0 




FIG. 5: Spectrums of a mobile hole for a 12 x 12 lattice in 
the BHM. The small dashed lines in (a) tag the dispersion of 
a free particle; in (d) they outline a dispersion of the form: 
E(p) = J|| (cos(p)x + cos(py)) 2 . 



Hence and as confirmed by Fig. [5] the dispersion rela- 
tion of a single hole in the BHM continuously deforms 
from the strong coupling form of Eq. i|31|) to that of a 
doped hole in a planar antiferromagnet (see Eq. I|32|)). 
Hence as a function of Jx/Ju there is a point where the 
effective mass (as defined by the inverse curvature of the 
dispersion relation) at p — (jr, it) diverges. Upon inspec- 
tion of the data (see Fig. [3J), the point of divergence of 
the effective mass is not related to the magnetic quan- 
tum phase transition and since it occurs slightly below 
• This crossover between a dispersion with min- 
imum at p — (-7T, 7r) and minimum at p = (tt/2, it/2) with 
a crossover point lying inside the AF ordered phase is 
also documented in Ref. |16|. 

The above argument can not be applied to the KNM, 
since the J± / Jn — point is macroscopically degenerate 
and hence is not a good starting point to understand the 
weak-coupling physics. Clearly the same holds for the 
KLM and UKLM. Inspection of the spectral data deep 
in the ordered phase of the KNM (see Fig. [U) shows that 
the maximum of the dispersion relation is still pinned at 
p = (-7T, 7r) such that the strong coupling features stem- 
ming from Kondo screening is still present at weak cou- 
plings. For the KNM and up to the lowest couplings we 
have considered the effective mass at p — (it, ir) increases 
as a function of decreasing coupling strength but does 
not seem to diverge at finite values of Jl/Jm. Precisely 
the same conclusion is reached in the framework of the 
KLM Q and UKLM 0. 



IV. QUASI PARTICLE RESIDUE 

In this section we turn out attention to the delicate issue 
of the quasiparticle residue in the vicinity of the magnetic 
quantum phase transition. We first address this question 
within the framework of the the mean-field model of Eq. 
(|29[) and compute the single particle Green's function 
within the framework of the self-consistent Born approx- 
imation. In a second step, we attempt to determine the 
quasiparticle residue directly from the Monte Carlo data. 



A. Analytical Approach 

Here we restrict our analysis to the BHM. and return 
to the Hamiltonian l|29|l . The coupling between the hole 
and magnons g(p,q) reads: 

g(p,q) = g a (p,q) + 9b(p,q)- 

We identify the two coupling constants with the processes 
that are shown in Fig. [7| g a (p,q) is proportional to the 
hopping matrix element and hence describes the coupling 
of a mobile hole to magnetic background, whereas gb(p, q) 

(2) 

is proportional to Jm and describes the coupling of a 
hole at rest with the magnons. Our calculations give the 
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FIG. 6: Spectrum of the KNM for a 12 x 12 lattice. The 
dashed lines tag the dispersion of a free particle. 



following momentum dependent coupling strengths: 
ts 

g a (p,q) = 7 =(-f{p + q)u(q)+j(p)v(q)) (33) 

V jV 



9bip,q) 



r f{q)(u(q)+v(q)) 



(34) 



where 7(g) = 2(cos((fc) + cos(q y )). We concentrate on 
the coupling to critical magnetic fluctuations and hence 
set q — Q and place ourselves in the proximity of the 
quantum phase transition, on the disordered side. In 
this case O(Q) — > and the coherence factors (see Eq. 
(12(it ) are both proportional to O(g) - 2 . Since furthermore 
lip + Q) = — 7(p) one arrives at the conclusion that 



.1 



i \ 

1 J 



f • ... ' • 

f f „ f f 



1 J 



1 J 



FIG. 7: Two possible processes where the hole can couple to 
magnons: (a) The hole moves to a next neighbor, (b) The 
hole is at rest. 



9a(p,Q) vanishes at the critical point. There is hence 
no coupling via process (a) to critical fluctuations. In 
other words process (a) couples only to short range spin 
fluctuations. On the other hand in the vicinity of the 
critical point gb scales as gbip,q) oc O(g) - 5 so that we 
can only retain this term to understand the coupling to 
critical fluctuations. Summarizing we set: 



gip,q) -> 9b(q) « 



1 



VO(g) 



(35) 



for the subsequent calculations. It is intriguing to note 

(2) 

that in this simple approximation gb(q) scales as Jy , 
which is strictly speaking null in the KNM. However, 
such a coupling should be dynamically generated via an 
RKKY-type interaction. 

With the above couple the first order self energy diagram 
for wave vectors satisfying e(p) = e(p + Q) shows a loga- 
rithmic divergence as a function of the spin gap. Hence 
we have to sum up all diagrams. We do so in the non- 
crossing or self-consistent Born approximation which in 
the T — limit boils down to the following set of self- 
consistent equations. 



E(p,w) 



l^ 3 2 (p,g)G(p 
g 

1 



g,w-0(g)) 



(36) 



Here we use a magnon dispersion relation of the form 
O(g) = VA 2 + v* (1 + 7(g)/4) with 7 (g) = 2(cos( fe ) + 
cos((?j,)), which agrees in the limit q — * Q = (tt, tt) with 
the form of Eq. (|27fl . Iterating the Green's function up 
to the 15th order to ensure convergence, we calculate the 
spectrum, pip, u) = ^Im[G(p, oS)\ via the imaginary part 
of the Green's function and compute the quasi-particle 
residue (QPR) at the first pole of the spectrum. 



Z(p) = 



1 



(37) 



Figure |H1 shows the QPR for p = (tt,tt) as a function of 
linear length L of the square lattice for different values of 
the spin gap A. The large-L limit is indicated by a line. 
Figure El plots the quasiparticle weight as a function of 
the spin gap for hole momenta p = (^, ^), (0, tt), (w, tt). 
For hole momenta satisfying e(p) = e(p + Q) ( p = (f , f ) 
and p — (0, tt) ) there is no energy denominator pro- 
hibiting the logarithmic divergence of the first order self- 
energy and the QPR shows an obvious decrease right up 
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FIG. 8: Self consistent Born approximation: QPR as a func- 
tion of linear lattice size, L, for different spin gap energies, 
A. 
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FIG. 9: Self consistent Born approximation: QPR in the 
vicinity of the quantum critical point for selected hole mo- 
menta (a) in a linear plot and (b) in a double logarithmic 
plot. A corresponds to the spin gap. 
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FIG. 10: Green's function in the vicinity of the phase transi- 
tion ( J±/J\\ = 2.5) for a static hole and various lattice sizes in 
the BHM (a) on a logarithmic plot and (b) on a plot where we 
adjusted the chemical potential in such a way that the Green's 
function converges to a constant value. Within the error bars 
and for lattice sizes greater then 12 x 12 there is no size scal- 
ing recognisable. (Inverse temperatures: /3J« — 30 (L — 12), 
PJ n = 50 {L = 16)), 0J n = 70 (L = 20); ArJ,, = 0.02) (c) 
QPR in the vicinity of the quantum critical point. 



to a complete vanishing at the critical point. Further- 
more, the data is consistent with 2 oc \fK. The case 
p = (tt, tt) is more complicated since e(p) ^ e(p + Q). In 
first order, the self-energy remains bounded. The scat- 
tering of the hole of Q — (tt, tt) magnons leads to the pro- 
gessive formation of shadow bands as the critical point 
is approached such that at the critical point, the relation 
EW(p) = EW(p + Q) holds. This back folding of the 
band can lead to the vanishing of the QPR when higher- 
order terms are included. Although the SCB results show 
a decrease of the QPR in the vicinity of the critical point, 
they are not accurate enough to answer the question of 
the vanishing of the QPR at this wave vector. 



B. QMC approach 

As shown in section [H] we can extract the QPR from the 
asymptotic behavior of the Green's function. We first 
concentrate on the static hole in the BHM for which the 
QMC data is of higher quality that for the dynamic hole. 
Fig. IlUb plots the Green's funtion as a function of lattice 
size at J±/J\\ = 2.5. As apparent within the consid- 
ered range of imaginary times no size and temperature 
effect is apparent. We fit the tail ( 5 < tJ» < 6) of the 
Green's function to the form Ze~ TfJ - and plot in Fig. IlUb 
G(T)e T/i . In the large imgaginary time limit this quan- 
titiy converges to the QPR Z. The so obtained value 
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FIG. 11: r)h (see Eq. 1381 ') as a function of J±/J\\ for a static 
hole in the BHM. 



of Z is plotted for values of J±/J\\ across the magnetic 
quantum phase transition. As apparent no sign of the 
vanishing of the QPR is apparent as we cross the quan- 
tum critcal point. 

Our QMC data allows a different interpretation. Follow- 
ing the work of Sachdev et al. 0] we fit the imaginary 
time Green function to the form: 



G(t) oc t Vh exp (— r/x) 



(38) 



in the the range 2.0 < tJh < 6.0 as done in Ref. |18j . 
Clearly, if rjh > then the QPR vanishes. Our results 
for rjh are plotted in Fig. ^] At J±/J\\ = 2.5 our result, 
i]h = 0.0875±0.0085 compares very well to that quoted in 
Ref. El, % = 0.087 ± 0.040. The fact that the result of 
Ref. [18j is obtained on a 64 x 64 lattice and ours on 20 x 
20 confirms that for the considered imaginary time range, 
size effects are absent. Given the above interpretation of 
the data, Fig. suggests that QPR of a static hole 
vanishes for all J±/J\\ < (J±/J\\) — 2.5. 
The choice of the fitting function reflects different order- 
ing of the limits t — * oo and N — > oo. On any finite size 
lattice the QPR is finite and hence it is appropriate to fit 
the tail of the Green's function to the form Z(N)e~ TfJ " to 
obtain a size dependend QPR, and subsequently take the 
thermodynamic limit. This strategy has been used suc- 
cessfully to show that the QPR of a doped mobile hole in 
a one dimensional Heisenberg chain vanishes . On the 
other hand, the choice of Eq. I|38fl for fitting the data im- 
plies that we first take the thermodynamic limit. Only in 
this limit, can the assymptotic form of the Green's func- 
tion follow Eq. Q38[) with % ^ 0. The fact that both pro- 
cedures yield different results sheds doubt on the small 
imaginary time range used to extract the quasiparticle 
residue. In particular, using data from tJ\\ = 2 onwards 
implies that we are looking at a frequency window around 
the lowest excitation of the order uj/J~ 0.5. Given this, 
it is hard to resolve the difference between a dense spec- 
trum and a well defined low-lying quasiparticle pole and 
a branch cut. 

We conclude this section by presenting data for a mobile 
hole in the BHM (see Fig. [T2J) and KNM (see Fig. IT3|) . 
Recall that in our simulations we restrict the motion of 




FIG. 12: Green's function of a dynamic hole (p = (vr,7r)) in 
the BHM for a 12 x 12 lattice at J±/J\\ = 2.4. (b) QPR in 
the vicinity of the quantum critical point. 



the hole to a single plane. The data for the QPR in the 
above mentioned figures stem from fitting the tail of the 
Green's function to the form Ze~ TM . The fit to the form 
of Eq. H38fl yields values of rjh which within the error 
bars are not distinguishable form zero. 



V. CONCLUSION 

We have analyzed single hole dynamics across magnetic 
order-disorder quantum phase transitions as realized in 
the Kondo Necklace and bilayer Heisenberg models. The 
hole motion is restricted to the upper layer as appro- 
priate for interpretation of the data in terms in Kondo 
physics. Both models have identical spin dynamics since 
the quantum phase transition is described by the 0(3) 
three-dimensional sigma model E3- On the other hand 
the single hole dynamics shows marked differences. In 
the strong coupling limit, deep in the disordered phase, 
the ground state of both models is well described by a di- 
rect product of singlets between the layers. This Kondo 
screening leads to a single hole dispersion relation with 
maximum at p = (tt, tt). In the Kondo Necklace model, 
where the spin degrees of freedom on the lower layer in- 
teract indirectly through polarization of spin on the up- 
per layer (RKKY type interaction), the single hole dis- 
persion preserves it's maximum at p — (tt,tt) down to 
arbitrarily low interplanar couplings. This situation is 
very similar to the Kondo Lattice model of Eq. QJ. In 
this case, down to J jt — 0.2, substantially below the 
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FIG. 13: (a) Green's function G p (t) in the Kondo Neck- 
lace model at p — (w,ir) for J±/J\\ = 1.0,1.5,2.0. (b) Ex- 
tracted values for the QPR. The critical point is localised at 
(J±/J\\)c = 1-360 ±0.017. 



magnetic phase transition J c /t = 1.45 ± 0.05, the max- 
imum of the the hole dispersion is pinned at p = (jr, 7r) 
and the effective mass at this p-point tracks the single 
ion Kondo temperature 0]. We note that this result 
is not supported by recent series expansions which show 
that there is a critical value of the coupling where the ef- 
fective mass diverges |2jj. Hence the interpretation that 
in both the Kondo necklace and Kondo lattice models, 
the localized spins remained partially screen down to ar- 
bitrarily low values of the interlayer coupling. In other 
words, signatures of strong coupling physics in the single 
hole dispersion relation is present down to arbitrary low 
interplanar couplings. 

In the bilayer Heisenberg model where there is an in- 



dependent energy scale coupling the spins on the lower 
layer, the situation differs. At values of J± < J± c the 
maximum of the single hole dispersion relation drifts to- 
wards p — (7r/2,7r/2) and the dispersion relation evolves 
continuously to that of a single hole doped in a planar 
antiferromagnetic |2l|. Hence the interpretation that at 
weak couplings, Kondo screening in this model is com- 
pletely suppressed. In other words, the small but fi- 
nite J± results can be well understood starting form the 
Jj_ = point. 

We have equally, analyzed the quasiparticle residue 
across the magnetic order-disorder transition. In the dis- 
ordered phase using a bond mean-field approximation, 
there are two processes in which the hole couples to mag- 
netic fluctuations (see Fig. [7J: i) The hole propagates 
from one lattice site to another thereby rearranging the 
spin background. In the proximity of the critical point, 
and still within the bond-mean field approximation those 
processes do not couple to long range Q = (-7T, 7t) mag- 
netic fluctuations. A very similar result is obtained in 
the ordered phase 01 • h) In bilayer models the hole can 
remain immobile and the spin in the lower layer can flip. 
Those processes couple to critical magnetic fluctuations. 
Within a self-consistent Born approximation, this drives 
the quasiparticle residue to zero both for a static hole and 
mobile hole with momenta p satisfying e(p + Q) = e(p). 
We have attempted to confirm this point of view with 
Monte Carlo simulations. Within our quantum Monte 
Carlo approach, where the accuracy of the single particle 
Green's function at large imaginary times is limited, we 
have found no convincing evidence of the vanishing of the 
quasiparticle residue both for a static and a mobile hole. 
Furhter work and algortihmic developments are required 
to clarify this delicate issue. 
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